Electron density distribution and screening in rippled graphene sheets 
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Single-layer graphene sheets are typically characterized by long-wavelength corrugations (ripples) 
which can be shown to be at the origin of rather strong potentials with both scalar and vector com- 
ponents. We present an extensive microscopic study, based on a self-consistent Kohn-Sham-Dirac 
density- functional method, of the carrier density distribution in the presence of these ripple-induced 
external fields. We find that spatial density fluctuations are essentially controlled by the scalar 
component, especially in nearly-neutral graphene sheets, and that in-plane atomic displacements 
are as important as out-of-plane ones. The latter fact is at the origin of a complicated spatial distri- 
bution of electron-hole puddles which has no evident correlation with the out-of-plane topographic 
corrugations. In the range of parameters we have explored, exchange and correlation contributions 
to the Kohn-Sham potential seem to play a minor role. 



PACS numbers: 71.15.Mb,71.10.-w,71.10.Ca,72.10.-d 

I. INTRODUCTION 

Graphene is a recently isolated material composed of 
carbon atoms arranged in a truly two-dimensional (2D) 
honeycomb lattice^™. States near the Fermi energy of a 
graphene sheet are described by a massless Dirac equa- 
tion which has chiral states in which the honeycomb- 
sublattice pseudospin is aligned either parallel to or op- 
posite to momentum. The Dirac-like wave equation and 
the existence of this spin-l/2-like quantum degree-of- 
freedom have a number of very intriguing implications 
on the properties of this material, most of which have 
been reviewed in the literature mentioned above. 

Graphene has been shown to possess a wealth of tanta- 
lizing electronic, mechanical, and optical properties and 
might well become the material that will replace silicon 
in the next generation devices^. Current exfoliated sam- 
ples however suffer from a limited mobility, with typi- 
cal values around 10.000 — 20.000 cm 2 /(Vs): the main 
source of disorder which is behind these numbers is not 
yet completely understood and represents a substantial 
obstacle against the quest for fundamental physical ef- 
fects and the development of functional devices. The 
mechanism which is limiting the mobility of the current 
(exfoliated) samples is actually one of the controversial 
topics in this field of research. Two "schools of thought" 
can be roughly identified: (i) one which ascribes the main 
limiting mechanism to charged impurities located in the 
(Si02) substrate™^, and (ii) one which instead relies on 
other scattering mechanisms, such as quenched rippled^, 
which are also long-range in nature. Ripples have been 
seen in suspended membrane d 14 * 15 * and also in flakes de- 
posited on substrates 16 19 and have been studied theoret- 
ically by Monte CarlcP^EU anc [ molecular dynamic d 22 * 23 * 
simulations. 

The controversy is enriched by several experiments 
which have targeted the role of disorder in exfoliated 



sampled^H3ll In particular, Bolotin et alW^ and Du et 
alW* have observed a drastic increase in mobility in sus- 
pended samples, in agreement with a scenario in which 
charged impurities in the substrate are the main source 
of scattering. On the other hand, Ponomarenko et alW* 
have studied exfoliated samples deposited on various sub- 
strates and found a rather weak dependence of the mobil- 
ity on the type of substrate. The authors of Ref. [29] have 
also studied transport in flakes embedded in media with 
very high dielectric constants, such as glycerol, ethanol, 
and water, and measured only a small increase in mo- 
bility. This experimental study seems thus to suggest 
that charged impurities are not necessarily the primary 
source of scattering in current samples. Whatever the 
leading sources of disorder are, it is of utmost importance 
to understand how well or poorly these are screened by 
electrons in graphene. 

The induced carrier density in graphene sheets sub- 
jected to the long-range potential of one or many charged 
impurities, in the absence or in the presence of electron- 
electron interactions, has been extensively studied theo- 
reticall}EM^ to the best of our knowledge, similar mi- 
croscopic studies in the presence of corrugations have not 
yet appeared. The aim of this article is to cover this 
gap: we present extensive self-consistent fully- quantum- 
mechanical calculations of the electronic density profiles 
of massless Dirac fermions in the external scalar and vec- 
tor potentials created by the corrugations. Our main 
findings can be summarized as follows: (i) the spatial 
density fluctuations induced by the ripples are almost 
entirely controlled by the scalar potential, especially in 
graphene sheets that are close to average neutrality; (ii) 
the contributions to the scalar and vector potentials due 
to in-plane atomic displacements are as large as those due 
to out-of-plane ones; and (iii) exchange and correlation 
contributions to the effective scalar (Kohn-Sham) poten- 
tial seem to play a minor role in determining the shape of 
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FIG. 1: (Color online) Three-dimensional plot of the corru- 
gated graphene sample used to calculate the average displace- 
ments shown in Fig. [2] and the scalar and vector potentials 
shown in Fig. [3] 



the ripple-induced electron- hole puddles, at least in the 
range of parameters we have analyzed. 

This manuscript is organized as follows. In Sect. [Il] we 
discuss in detail how we have calculated scalar and vector 
potentials starting from a corrugated graphene sheet. In 



Sect. Ill we introduce the theory and the numerical pro- 
cedure we have used to calculate the induced carrier den- 
sity in the presence of the ripple-induced potentials and 
present our main numerical results. Finally in Sect. IV 
we draw our main conclusions. Appendix [A] reports some 
technical remarks on the calculation of the density in- 
duced by a purely vector potential within linear-response 
theory. 



II. FROM RIPPLES TO SCALAR AND 
VECTOR POTENTIALS 

The aim of this Section is to describe how we have com- 
puted the scalar and vector potentials associated with 
ripples. For definiteness we focus our attention on rip- 
ples generated by thermal fluctuations^^. The proce- 
dure we have followed is however completely general and 
applies to any type of ripples, independently of the mi- 
croscopic, intrinsic or extrinsic, mechanisms that lie at 
their origin. 



A. Microscopic calculation of the average 
displacements 

In what follows we consider a specific realization of a 
corrugated graphene sheet at a temperature T = 300 K, 
computed with a Monte Carlo simulation as in Ref. [21] 
In Fig. [I] we show the three-dimensional bond structure 



of the sample, which contains 19504 atoms and fulfills 
periodic boundary conditions in the simulation box. 

The computation of the corrugation-induced scalar and 
vector potentials that we will carry out in Sect. |IIB] be- 
low requires the knowledge of the displacements {ui} of 
the atomic positions {r-} in the sample (i is the atomic 
label) with respect to a flat reference distribution {r^}. 
The latter is defined by applying a dilation/contraction 
to the honeycomb lattice at T = 0. More precisely, we 
first make sure that the positions, tcm and r' CM , of the 
center-of-mass of the two distributions coincide, and use 
in the following the displaced vectors r — » r — vcm- We 
then dilate/contract the honeycomb lattice at T = to 
compensate for the variation of the carbon-carbon bond 
length produced by the finite temperature. The coef- 
ficient A in the transformation r — » Xr is obtained by 
averaging the ratio A^ = |^|/|r^| over all the atoms i 
such that \ri\ > 50.0 A. The latter restriction reduces 
the impact of the fluctuations of the atomic positions, 
produced by the ripples, but does not affect the com- 
putation of the overall stretch/compression produced by 
the temperature. We find A c± 0.998 (< 1: the effect of 
temperature in this range is indeed to reduce the carbon- 
carbon bond length^). The variance of {A^} is of order 
10 -3 , hence the stretch induced by the temperature is the 
dominant contribution of the atomic displacements from 
the positions of the bare honeycomb lattice. In other 
words, to prepare a sensible reference distribution it is 
essential to perform the aforementioned stretch, even if 
the factor A is close to unity. 

Finally, we make sure that the sample and the reference 
distribution are not globally rotated with respect to each 
other. We compute the average angular displacement 
vector 







Vi X Ti 



(1) 



where the summation is restricted to the atoms such 
that the cosine of the angle between r[ and is larger 
than 0.9. In the analyzed sample the modulus of <fi is 
of order 10 -3 , hence we conclude that the sample and 
the reference distribution are properly aligned. We are 
now in position to compute the displacement vectors 
Ui = r\ — Ti'. thanks to the above mentioned prepara- 
tion procedures, these will be free of artificial systematic 
trends and will provide us with an accurate local descrip- 
tion of the ripples. 

As we solve for the electronic density on a square mesh 
in the simulation box (see the description of the method 
in Sec. IIIB), the knowledge of the displacement of each 
atom is superabundant. For this reason we average the 
atomic displacements over square patches defined on a 
square mesh. To show that this averaging yields indeed a 
correct modeling of the physical system, we observe that 
the problem possesses three length scales: (i) graphene's 
lattice constant a = aoV^ ~ 0.25 nm (here ao = 1.42 A is 
the carbon-carbon distance), (ii) the length scale A s of the 
spatial structures in the specific sample shown in Fig. [I] 
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FIG. 2: (Color online) Average displacements u(r) calcu- 
lated as discussed in Sect. |II A| The color scale represents 
the z component of the average displacements, varying from 
—3.0 A (blue) to +3.0 A (red). The arrows, whose length has 
been multiplied by a factor ten for better visibility, represent 
the in-plane components of the average displacements. 



which is of the order of several nanometers (A s « 8 nm); 
and (iii) the spatial resolution A res which we have in our 
continuum-model electronic structure calculations [see 



Eq. (28)]. For a sample of roughly 22 nm x 22 nm, as 
the one shown in Fi g, [I] A res « 1.5 nm (see discussion 



IIIB[ ). Since A s ^> A res ^> a, the struc- 



below in Sect. 

tures in Fig. [l] are properly resolved by the mean dis- 
placement vectors u(r), obtained by averaging the micro- 
scopic displacements over square patches of area « A 2 es . 
The result of this averaging procedure for the sample in 
Fig. [I] is shown in Fig. [5] where we have plotted u(r) as 
calculated on a square mesh with 32 x 32 points. We 
remark that the in-plane displacements undergo strong 
variations between neighboring patches as a consequence 
of the fact that even the in-plane displacements of neigh- 
boring atoms in the sample do not present signatures of 
local correlations. 

We now proceed to discuss how we have calculated the 
deformation tensor and the corrugation-induced scalar 
and vector potentials. 



where Uij (with 
tensor, 



5 j £ { x ->y}) is the usual deformation 



dm du 



dXj 
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du k du k 
dxj dxn 



(4) 



Here Ui = Ui(r) with i G {x, y, z} are the Cartesian com- 
ponents of the average displacements. For the coupling 
constant g\ we have used two v alues, g\ = 3 eV and 
gi = 16 eV (the latter valu d 45 * 4 ^ , which is based on old 
transport data on graphite sample, seems largely overes- 
timated 4 ^), while 
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where f3 = — dlog (70) /Slog (ao) ~ 2, 70 « 2.7 eV is the 
nearest-neighbour hopping parameter, and 



1 Ms 
y/2 B 



(6) 



For the shear fi s and bulk B moduli we have used the 
recently calculated valued, /i s = 9.95 eV A -2 and B = 
12.52 eV A -2 , at a temperature T = 300 K. We thus 
find that k, ~ 0.56 at this temperature. 

In Fig. [3] we illustrate scalar and vector potentials cal- 
culated using Eqs. above. While performing the 
calculation of V± and V2 we have noticed that the deriva- 
tives of the average in-plane displacements u± are of 
0(1O -2 ), while the derivatives of the out-of-plane dis- 
placements u z are much bigger, O(10 _1 ). However, in 
the deformation tensor Q the latter enter only quadrat- 
ically. We thus conclude that the contributions from 
in-plane and out-of-plane displacements are both of the 
same order, O(10 _1 ). As a result, no evident correla- 
tions link the out-of-plane topographic corrugations [i.e. 
the distribution of the out-of-plane average displacements 
u z (r) shown in the color map in Fig. E with the scalar 
and vector potentials illustrated in Fig. 



III. KOHN-SHAM-DIRAC 
DENSITY-FUNCTIONAL CALCULATIONS 



B. The deformation tensor and the 
corrugation-induced scalar and vector potentials 

We have calculated scalar V\ and vector V2 = A x — 
iA y potentials according to the standard formulas of the 
theory of elasticity* 45 1 46 1 : 



and 



Vl = gi(U XX + Uyy) 



V 2 = g2(u xx - Uyy + 2iu xy ) , 



(2) 



(3) 



In this Section we present an approximate self- 
consistent microscopic theory for the carrier density dis- 
tribution in the corrugation-induced scalar and vector 
potentials shown in Fig. [3j 



A. Approximate Kohn-Sham-Dirac theory for 
corrugated graphene sheets 

We have generalized the Kohn-Sham-Dirac (KSD) the- 
ory described in Ref. [40] to deal with situations in which 
the massless Dirac fermion liquid is subjected to a space- 
dependent vector potential A(r) (the vector potential 
introduced below has the physical dimensions of energy) 
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FIG. 3: (Color online) Left panel: color plot of the scalar potential Vi(r) (in units of meV) calculated using Eq. (2} with 
gi = 3 eV. Central panel: the real part of the potential V^ir) (in units of meV) calculated using Eq. Right panel: the 
imaginary part of the potential V2(r) (in units of meV). 



which changes smoothly over many lattice constants. In 
this limit the induced density 5n(r) can be calculated by 
solving the following single-spin single-valley KSD equa- 
tion: 

{a ■ [vp + A(r)] + t a V KS (r)} $ A (r) = e\$\(r) . (7) 

Here a is a 2D vector constructed with the 2x2 
Pauli matrices g\ and a 2 acting in pseudospin space, 
v = 3joclq/(2H) « 10 6 m/s is the bare Fermi velocity, 
p = —ihV ri t a is the 2x2 identity matrix in pseudospin 
space, and the Kohn-Sham potential, 



V KS (r) = Vkt(r) + AVn(r) + V xc (r) , 



(8) 



is the sum of the external scalar potential V^ xt (r), the 
Hartree potential, and the scalar exchange-correlation 
potential. For A = Eq. (7k reduces to the KSD equa- 
tion introduced in Ref. |40J Note that Eq. ^ neglects 
exchange-correlation corrections to the vector potential 
A, which are beyond the scope of the present paper and 
which will be addressed in a subsequent publication. 

The ground-state density n(r) is obtained as a sum 
over the KSD spinors &\(r): 



n(r)= 9 ^[ A \r)\ 2 + \ V >[ B \r)\ 2 ]f(e x ) 



(9) 



where the factor g = g s g v = 4 is due to valley and 
spin degeneracies, {ip^\r),<j = A, B} are the pseu- 
dospin (sublattice) components of the spinor &\(r), and 
f(x) = {exp [(x — ii)/{k&T)\ + is the usual Fermi- 
Dirac thermal factor at a chemical potential \i = fi(T). 
Equation ([9| is a self-consistent closure relationship for 
the KSD equation (|7|), since the Kohn-Sham potential 
Vks(v) is a functional of the ground-state density n(r). 

In the absence of any source of external scalar and mag- 
netic fields, the scalar V ext (r) and vector A(r) potentials 
are solely determined by the corrugations: 



f V ext (r) = V 1 (r) 
{ A(r) = (5Re V 2 (r), -9m V 2 {r)) 



The Hartree potential is given by 



AVk(r) = / dV e " 5n(r') 
J e\r-r'\ 



(10) 



(11) 



where e is an average dielectric constant 



ei + e 2 



(12) 



Here e\ and €2 are the dielectric constants of the media 
above and below the graphene flake. For example e ~ 2.5 
for graphene placed on Si02 with the other side being 
exposed to air, while e « 1 for suspended graphene. The 
quantity Sn(r) = n(r) — no is the local density measured 
relative to a "background" value, no, which is defined by 



n 



2 

A 



(13) 



Here 2/^4o is the density of a neutral graphene sheet, 
^4o = SV^o/^ ~ 0.052 nm 2 being the area of the unit 
cell in the honeycomb lattice, and n c is the spatially av- 
eraged carrier density, which can be positive or negative 
and controlled by gate voltages. 

The third term in VKs(r), V xc (r), is the scalar 
exchange-correlation potential. This is a functional of 
the ground-state density, which is known only approxi- 
mately. Following Ref. [40] we employ the local-density 
approximation (LDA), 



Kc(r) 



= V hom (7 



01 , 
' I n— »n c (r 



(r) 



(14) 



where v^° m (n) is the T — exchange-correlation poten- 
tial of a uniform 2D liquid of massless Dirac fermiond 40 * 50 * 
with carrier density n. v^ m (n) * s related to the ground- 
state energy per excess carrier 5e xc (n) by 



d[nfe xc (n)] 
dn 



(15) 



The carrier density n c (r) is the density relative to that 
of a uniform neutral graphene sheet: 

2 

n c (r) = n(r) — = n c + Sn(r) . (16) 

Ao 

The expression used for Se xc (n) depends on the zero-of- 
energy, which is normally^ chosen so that v^° m (n — 0) = 
0. 



5 



B. Technical remarks on the method of solution 

In order to solve Eq. ^ we have followed the same 
technique adopted in Ref. 40, i.e. we use a square sim- 
ulation box of size L x L with periodic boundary con- 
ditions and conveniently expand the spinors $>\(r) in a 
plane- wave basis. We discretize real space restricting r 
to a square mesh = (i5, j5), with z, j — 1, . . . , N. Here 
S = L/N is the spacing of the mesh. Fourier transforms 
f(k) of real-space functions f(r) are calculated by means 
of a standard fast-Fourier-transform algorithm 51 that al- 
lows us to compute / on the set of discrete wavevectors 
h 



k 



2tt 



(17) 



with —N/2 < n Xj i,n y j < N/2 (or, equivalently, < 
n x ,i,n yJ < N). 

In momentum space Eq. ([7| reads 

^(k\{a•[vp^A(r)}^l a V K s(r)}\k , }^ x (k , ) = e x $\(k) , 

k' 

(18) 

and the problem is thus mapped into the diagonaliza- 
tion of the KSD matrix = (k\{a • [vp + A(r)] + 

I^^sM }!&')• The matrix elements in Eq. ( p~8| ) can be 
computed either analytically or numerically. More specif- 
ically, the matrix elements of the kinetic Hamiltonian are 
given by 



(k\ va • p \k'} = hvcr • k' 5k^> 



(19) 



The matrix elements of the Hartree term are given by 



<fc|AVk(r)|fc') 



27re 2 



elfe — fc 



- Sn(k - k') , 



(20) 



where Sn(k) = n(k) — rio^fc,o is the Fourier transform of 
the charge neutral density 5n(r), introduced above. 

The matrix elements of the external, vector, and 
exchange-correlation potentials can be calculated numer- 
ically from 



<fc|/(r)|fc'> 



1 

L2 



d 2 r f(r) 



-i{k — k')-r 



(21) 



where f(r) is either V ext (r), V xc (r), A x (r), or A y (r). 

In practice the diagonalization of the KSD matrix 
%k S jP requires the introduction of a momentum space 
cut-ofT 40 , k Xj i,k y j G [— fc c , + fc c ], which does not exceed 
the Brillouin-zone boundary defined by our real-space 
discretization: k c < tt/5. k c defines the range of mo- 
menta used in the expansion of the Hamiltonian H^|P 
and thus defines its dimension du' 



2 x 2 x 



Lk c 
~2~n 



(22) 



The factor of 2 here is due to the sublattice pseudospin 
degree-of-freedom. Given a value of k c the Kohn-Sham- 
Dirac matrix T-L^^P has dn eigenvalues, labeled by the 
discrete index A = 1 , . . . , dn . 

Let us consider a neutral-on- average graphene sheet 
(n c = 0) with areal extension L x L. The total number 
of electrons in such sheet is 



Ni 



real 



z r 2 
~T~ X L 



(23) 



The total number of electronic states available in our cal- 
culations is gdn. To simulate a neutral-on- average sheet 
we clearly need half of these states: 



N 



simul 



^ x gdn 



, Lk r 
5 x(2x^ + l 



(24) 



In Ref. 00] the authors enforced the following condition 



N simul — , 



real 



(25) 



which physically means that all the electrons in the 7r- 
band are simulated. This leads to the relation 2L 2 / A® = 
g [2Lk c /(27r) + l] 2 which links the momentum-space cut- 
off k c and the size of the system L. This relationship is 
however too restrictive since one would need very large 
values of k c (much larger than those prescribed by the 
computational limit) to simulate flakes with an areal ex- 
tension of experimental interest^. Therefore, the re- 
quirement p5| ) severely affects the possibility of perform- 
ing quantitative predictions for large systems. There are 
also more physical reasons for lifting the requirement 
(25): the massless Dirac fermion modeP does not de- 
scribe all electrons in the 7r-bands but only a fraction 
rj <C 1 of them. We thus have decided to relax the con- 



straint (25) allowing A^ S i mu i ^ iV rea i, i.e. 



N 



simul 



real 



(26) 



with < rj <C 1. Letting rj be different from unity we 
can choose L and k c independently. The factor rj' can be 
tuned in order to fulfill Eqs. fl23|, plf, and fl26|: 



V 



4 L 2 



= g [2Lk c /(2ir) + if 



2L 2 



(27) 



For example, we can choose L « 22 nm (as in the case of 
Fig.[T]) and fix k c according to our numerical capabilities, 
say k c = 15 x (2n/L). Substituting these values for L and 
k c in Eq. (27), one obtains that the fraction of simulated 
electrons in this case is rj « 0.2, i.e. 20% of the electrons 
in graphene's 7r-band. We remark that the existence of 
a momentum space cut-off k c implies a minimum spatial 
resolution, 



2tt 

k r 



(28) 



which in this case would be A res ~ 1.5 nm, and thus suffi- 
cient to resolve rather short-wavelength spatial structures 
in the induced carrier density. 
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The arguments above can be readily generalized to the 
case of a doped graphene sheet (n c ^ 0): in this case 
Eq. (f26l) reads 



simul 



^+n c L 2 = r/2^+n c L 2 . (29) 



We clearly see that even at finite doping we can arbitrar- 
ily choose L and & c , with a fraction of simulated electrons 
which is still given by Eq. (27). 

Before concluding this Section we recall that the ex- 
change and correlation potential v^° m (n) introduced in 
Sect. |IIIB] depends on carrier density n c through the di- 
mensionless quantit}^ A = /c m ax/^F, where /c max is an 
ultraviolet cut-off and /cp = \/^\n c \/ g is the Fermi wave 
number. We take k max to be such that 



Sir 2 

gA 



(30) 



where rj is a dimensionless number, < rj < 1, which 
should be assigned a value according to the wave vec- 
tor range over which the continuum model describes 
graphene^. Thus, making use of Eqs. (27) and (30), 
we find 



A 



27/ 



A \n c 




2\n r \L 2 



(31) 



However, it is physically reasonable to identify r/ and rj' 
since they both refer, directly or indirectly, to the range 
of applicability of the massless Dirac fermion model to de- 
scribe electrons in graphene. Consequently, we see that, 
taking r/ = r/, A is independent of the choice of r/ while 
it depends on n c L 2 , i.e. on the average carrier density 
in units of 1/L 2 , and on the dimension of the KSD 
Hamiltonian (or equivalent ly on k c ). 
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FIG. 4: (Color online) Top panel: fully self-consistent elec- 
tronic density profile 5n(r) (in units of 10 12 cm -2 ) in a 
corrugated graphene sheet. The data reported in this fig- 
ure have been obtained by setting gi = 3 eV, a ee = 0.9 
(this value of a ee is the commonly used value for a graphene 
sheet on a Si02 substrate), and an average carrier density 
n c c± 0.8 x 10 12 cm -2 . Bottom panel: same as in the top 
panel but for a ee = 2.2 (this value of a ee corresponds to sus- 
pended graphene). 



C. Numerical results 



In Fig. [4p0] we report our main numerical results ob- 
tained from the self-consistent solution of the KSD equa- 
tion ([7]) with a momentum-space cut-off k c = 15 x (2tt/L). 
The induced density profiles depend on the strength of 
electron-electron interactions which is measured by the 
dimensionless fine-structure constant 



ehv 



(32) 



In Fig.|4]we illustrate the fully self-consistent electronic 
density profile Sn(r) in the ripple-induced scalar and vec- 
tor potentials shown in Fig. [3] By "fully self-consistent" 
we mean that Sn(r) has been obtained with the inclusion 
of both Hartree and scalar LDA exchange-correlation po- 
tentials. In this figure we have reported results for two 
values of graphene's fine structure constant, a ee = 0.9 
(graphene on Si02) and 2.2 (suspended graphene). We 



clearly see electron-hole puddles with a typical size of a 
few nanometers. 

In Fig. [5] we show one-dimensional cuts of 5n(r) for the 
same system parameters as in Fig.[4]to better address the 
separate role of Hartree and exchange-correlation poten- 
tials. From the top panel in Fig. [5] we clearly see what 
is the role of electron-electron interactions and screen- 
ing: the amplitude of the density fluctuations is indeed 
completely controlled by interactions. From the bottom 
panel we see how, for this particular set of parameters, 
scalar LDA exchange and correlations effects seem to be 
playing only a minor (quantitative) role. 

As in Ref. 40, it is interesting to compare the reduc- 
tion in the amplitude of density fluctuations seen in the 
top panel of Fig. [5] with what would be expected in a lin- 
ear screening approximation. Assuming that the biggest 
role is played by the scalar poten tial Vj (this assump- 
tion will be justified below in Sect. IIIC), within linear- 



response theory (LRT) the induced density change (in 
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FIG. 5: (Color online) Top panel: a one- dimensional plot of 
Sn(r) (as a function of x in nm for y — 11.3 nm) for the 
same set of parameters as in the lower panel of Fig. ^] Here 
we have reported data for noninteracting electrons (circles), 
data obtained including only the Hartree term in Eq. (pi 
(triangles), and data obtained including both Hartree and 
exchange-correlation potentials (squares) . Note that electron- 
electron interactions completely control the magnitude of den- 
sity fluctuations and that, on this scale, the data obtained 
including exchange-correlation effects (squares) are indistin- 
guishable from the data obtained with the inclusion of the 
Hartree potential only (triangles). Bottom panel: same as 
in the top panel but with the exclusion of data for noninter- 
acting electrons. Differences between data labeled by squares 
and by triangles can be seen on this scale. These differences 
are however only quantitative and not qualitative. 



Fourier transform) is given by 



(33) 



where Xo(q) is the static T = Lindhard function of a 
homogeneous noninteracting massless Dirac fermion fluid 
(see for example Ref . l50|) , 

/ \ / x 9Q „ ( 2k F \ gk F ( 2k F 

(34) 

and e(q) = 1 — v q xo(q) is the static random-phase- 




FIG. 6: (Color online) A one-dimensional plot of Sn(r) (as 
a function of x in nm for y — 11.3 nm) for the same set of 
parameters as in Fig. [5] Here we compare results based on the 
solution of Eq. |7| with electron-electron interactions treated 
at the Hartree level (triangles) with those based on linear- 
response theory (hexagons), Eqs. (33)- (36). Linear screening 



seems to describe very well the data. 



approximation dielectric function: 

/ \ , , Qtf . 7T f2k F \ q^ F (2k F 

e(q) = 1 + +g-a ee F G 

q 8 \ q J 2q \ q 



(35) 



Here ^(^f) = gk F /(27rhv) is the density-of-states at the 
Fermi level, v q = 2ne 2 /{eq) in the Fourier transform of 
the electron-electron interaction, ^tf = g&eek F is the 
Thomas-Fermi screening vector, and, finally, 



1 arcsin 

7T 



F(x) 

k G{x) = y/l-x 2 6(1 -x) 



(36) 



Note that F(x) = G(x) = for x > 1 (i.e. q < 2k F ). 
In Fig. [6] we show a comparison between the prediction 
of LRT, based on the Fourier transform of Eq. ( |33| ), and 
the non-linear screening result based on the solution of 
Eq. ([7]) with the Hartree potential only. We thus see 
that, maybe surprisingly, LRT explains the data quanti- 
tatively. 

In Fig. [7] we show fully self-consistent electronic density 
profiles obtained for a much larger value of the scalar 
g\ constant. These results have to be compared with 
those reported in Fig. [4] As expected, in the case g\ — 
16 eV the amplitude of the density fluctuations is much 
larger. A direct comparison has been reported in the 
one-dimensional cuts in Fig. [8j 

The dependence of the self-consistent density profiles 
on the doping level n c is shown in Fig. |9] from this plot, 
and especially from the inset, we see that the ampli- 
tude of the density fluctuations see m to s aturate slowly 
with increasing n c , as already founcPSEH in the case of 
self-consistent screening calculations in the presence of 
randomly-distributed charged impurities. 
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FIG. 7: (Color online) Same as in Fig.[4]but for gi = 16 eV. 
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FIG. 8: (Color online) A one-dimensional plot of the fully self- 
consistent 5n(r) (as a function of x in nm for y = 15.8 nm) ob- 
tained using gi = 3 eV (circles) or g\ — 16 eV (triangles). The 
other parameters are a ee = 2.2 and n c c± 0.82 x 10 12 cm -2 . 



Before concluding this Section we stress again that 
there is no evident correlation between the out-of- 
plane topographic corrugations and the spatial structures 
(electron-hole puddles) in the density profiles, as already 
pointed out in Sect. |IIB| This is highlighted in Fig. [l0| 



FIG. 9: (Color online) One-dimensional plots of the self- 
consistent density profiles (as functions of x in nm for y = 
21.1 nm) for different values of doping: n c — 0.8 x 10 12 cm -2 
(circles), n c — 3.96 x 10 12 cm -2 (triangles), and n c — 
3.17 x 10 13 cm -2 (squares). The data reported in this fig- 
ure have been obtained by setting g± = 3 eV and a ee = 2.2. 
The inset shows 5n(r) (in units of 10 12 cm -2 ) at a given point 
r in space as a function of the average carrier density n c (in 
units of 10 12 cm" 2 ). 



D. 



Self-consistent electronic density in the 
presence of a model ripple 



As emphasized in Sects. II B and III C[ in-plane and 
out-of-plane displacements have the same impact on the 
corrugation- induced scalar and vector potentials: this re- 
sults into complicated spatial patterns of the carrier den- 
sity with no immediate link with the topographic corru- 
gations. In this Section we present the self-consistent 
electronic density profile in the presence of a simple 
model ripple which exhibits displacements only in the 
z direction. 

For concreteness, following Ref. [54j we consider the 
following Gaussian out-of-plane displacement: 



u z {r) = A exp 



°rel 



b 2 



(37) 



where x re \ = x — L/2 and y re \ = y — L/2. The scalar and 
vector potentials can be easily computed from Eqs. Q 
and ([3]), leading to the following expressions: 



Vi(r) 



A 2 



2 

rel 



2/rel) eX P -2 



°rel 



Vrel 



and 



V2(r) 



(x 



hA v-rel + Wrel) exp . _ ^ 

'(39) 

The fully self-consistent density profile 5n(r) calcu- 
lated with the use of the potentials (38) and ( |39| ) is re- 
ported in Fig. [TT] These data show that when in-plane 
displacements are neglected the correlation between the 



b 2 



(38) 
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FIG. 10: (Color online) Three-dimensional plot of the fully 
self-consistent continuum- mo del Dirac-Kohn-Sham density 
profile reported directly on the corrugated graphene sample 
shown in Fig.[l] More precisely, the color-coding of the hexag- 
onal bonds labels the local value of Sn(r) shown in the two- 
dimensional color plot reported in the bottom panel of Fig. [4] 
Note that there is no simple correspondence between the out- 
of-plane topographic corrugations and the density profile. 
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FIG. 11: (Color online) Top left panel: color plot of the an- 
alytical scalar potential Vi(r) (in units of meV) reported in 
Eq. (38). The parameters used are: gi = 3 eV, A = 0.05L = 
1.1 nm, and b = 0.2L = 4.5 nm. Top right panel: real part of 
the analytical potential \^(r) (in units of meV) in Eq. (39). 
Bottom left panel: imaginary part of the potential V^(r) (in 
units of meV) in Eq. (39). Bottom right panel: fully self- 



consistent electronic density profile (in units of 10 12 cm -2 ) 
calculated in the presence of the scalar and vector potentials 
shown in the other panels. This numerical calculation has 
been performed using a ee = 2.2 and n c — 3.96 x 10 12 cm -2 . 



density profile and the topography of the corrugated 
graphene sheet [Eq. ( |37[ )] is much more transparent. Note 
that the oscillations in 5n(r) stem from the fact that the 
quantity \\7u z (r)\ 2 , which controls the scalar potential 
Vi, is maximal at \r\ « b. 



E. Comments on the density response to a purely 
vector potential 

A natural question might arise at this point: what is 
the relative role of V\ and V2 in determining the induced 
density 6n(r)7 In this Section we study the density re- 
sponse of a system of massless Dirac fermions to a purely 
vector potential. 

Let us begin for simplicity from a noninteracting sys- 
tem: in this case we can prove that Sn(r) = 0, inde- 
pendently of doping. This can be easily seen within the 
framework of LRT: in this case 



Sn(q) = Xn 3 i(o)Ai(q) 



(40) 



where Sn(q) and Ai(q) are the Fourier transforms of 
Sn(r) and A^r), and Xnj^o) = lim ^oXn/fe^) is a 
static linear-response function. It turns out (see Ap- 
pendix [A] for a formal proof) that 



(41) 



where Xnn(Q^) is the density-density response function 
of a noninteracting system of massless Dirac fermions 
(see for example Ref.l50land references therein). Because 
Xnn(Q,u) is well behaved in the static limit we immedi- 
ately find that Xnf (o) — 0. 

An identical conc lusion can be reached by invoking 
Furry' s theore m 55 * 56 ! , which applies independently of the 
strength of the external vector potential A (and thus also 
beyond the regime of applicability of LRT) and in the 
presence of electron-electron interactions. The theorem, 
however, is valid only for systems with an electron-hole- 
symmetric spectrum. We thus expect Sn(r) = only in 
the case of a neutral-on-average system, while we expect 
a finite induced density for a finite value of n c . 

We have checked these expectations numerically. We 
have performed calculations in the presence of the scalar 
V\ component only and compared the calculated induced 
density, 6ns (r), with that obtained in the presence of 
both scalar and vector potentials, (^totM- In Fig. 12 
we report the results for a ee = 0: we clearly see, espe- 
cially from the bottom panel, that even at finite average 
carrier density the amplitude of the spatial fluctuations 
induced by the vector potential only is rather small. 

Differences between 6ns (r) and ^nxoT(^) have been 
quantified by the value of the following dimensionless pa- 
rameter, 



\/\\6n TOT (r) - 6n s (r)[\ 
Vll^TO T (r)|| + Vll^s(r)|| ' 



(42) 
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FIG. 12: (Color online) Top panel: a one-dimensional plot of 
the noninteracting (a ee = 0) density profile Sn(r) (as a func- 
tion of x in nm for y = 12.3 nm) obtained solving the Dirac 
equation in the presence of both scalar and vector potentials 
(circles) or of the scalar potential only (triangles). Bottom 
panel: a one-dimensional plot of the noninteracting density 
profile 6n(r) (as a function of x in nm for y — 12.3 nm) ob- 
tained solving the Dirac equation in the presence of the vector 
potential only. The data reported here refer to g\ = 3 eV and 
n c — 3.96 x 10 12 cm -2 . From both panels we conclude that 
density fluctuations are largely controlled by the scalar po- 
tential. 



F. Electronic density in the presence of both 
ripples and charged impurities 

Before concluding we would like to briefly illustrate 
how the presence of the ripples modifies qualitatively the 
density landscape indu ced by a random distribution of 
charged impuritie d 4Q * 41 l In this Section we report nu- 
merical results based on the self-consistent solution of 
Eq. ([7]) in the presence of a scalar potential V GX t(r) given 
by: 



V ext (r)=V 1 (r) + V imp (r) . 



(44) 



Here Vj mp (r) is a scalar potential due to charged impu- 



ritie; 



40 



^mp(r) 



AW 



Ze 2 



d 2 



(45) 



where Ri are random positions in the supercell and d is 
the distance between the graphene sheet and the plane 
where the impurities are located. For simplicity, all 
charges have been taken to have the same Z in Eq. (45). 



In Fig. [13] we show fully self-consistent density pro- 
files of massless Dirac fermions subjected to the scalar 
potential of A/"i mp = 5 charged impurities: in the top 
panel we show Sn(r) calculated in the absence of rip- 
ples (gi =#2=0), while in the bottom panel we have 
included them. We clearly see how the smooth land- 
scape of electron-hole puddles in the presence of charged 
impurities only (top panel) is dramatically affected by 
the presence of corrugations (bottom panel), which in- 
duce additional spatial variations with a much smaller 
length scale (probably well below the current spatial ex- 
perimental resolution of probes like SET^ or STiVpl). 
Once again, we would like to emphasize that these small- 
wavelength carrier-density oscillations are due to a com- 
plicated interference between the effects of out-of-plane 
and in-plane atomic displacements. 



IV. CONCLUSIONS 



where 



||0(r)|| 2 = J d 2 r\0(r)\ 2 (43) 



is the usual L 2 norm. In the case n c = we find 
£ ~ 3 x 10 -4 , which is below our numerical precision 
(0.005): within the accuracy of the calculation thus 
^tot(^) = 5ns(r). In the calculations with finite car- 
rier density, however, we find much higher values of e: 
for n c ~ 3.96 x 10 12 cm -2 we find e ~ 0.02, while for 
n c ~ 3.17 x 10 13 cm -2 we find e ~ 0.03. 



In summary, we have presented quantitative calcula- 
tions of scalar and vector potentials induced by corruga- 
tions in single-layer graphene sheets. We have found that 
the contributions from in-plane and out-of-plane atomic 
displacements are both of the same order and that this 
does not lead to evident correlations between the out-of- 
plane topographic corrugations and the induced scalar 
and vector potentials. 

We have then used these potentials to calculate self- 
consistently the induced electronic density distribution 
in the presence of electron-electron interactions. To this 
end we have generalized the Kohn-Sham-Dirac theory of 
Ref. [40] to treat situations with spatial-dependent vec- 
tor potentials. We have discovered that spatial density 
fluctuations are largely controlled by the scalar poten- 
tial, especially in nearly-neutral graphene sheets, and 
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FIG. 13: (Color online) Top panel: fully self-consistent elec- 
tronic density profile (in units of 10 12 cm -2 ) calculated from 
the solution of Eq. |7| in the presence of A^i mp = 5 charged 
impurities with charge Z = +1 (donors). The white circles 
label the position of the charges on a plane located at a dis- 
tance d c± 2 nm from the graphene sheet. Bottom panel: same 
as in the top panel but in the presence of ripples too. The 
data reported here have been obtained by setting g\ — 3 eV, 
a ee = 0.9, and n c ~ 3.96 x 10 12 cm" 2 . 



that this creates complicated short-wavelength (a few 
nm) electron-hole puddles which do not exhibit evident 
correlations with the topography of the sheet. 

In the future we would like to investigate more deeply 
the role of the exchange-correlation corrections to the 
vector potential, especially in view of the fact that 
the exchange-correlation contribution to the scalar Kohn- 
Sham potential, which has been studied here, has been 
found to play a minor role. 
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Appendix A: Density response to a vector potential 

In this Appendix we demonstrate that within LRT an 
external vector potential does not induce density mod- 
ulations in a system of noninteracting massless Dirac 
fermions (MDFs). 

We consider the following Hamiltonian (h = 1 in this 
Appendix): 



H = Ho + H' 



(Al) 



where 



Ho = -iv J2 I d2r tkr)**? ■ V ^( r ) ( A2 ) 

a, (3=1 ^ 



is the MDF kinetic Hamiltonian and 

H'= [ d 2 r A(r,t) ■ j(r) , 



(A3) 



A(r,i) being a weak perturbing vector potential acting 
on the system. Here we have introduced the well-known 
MDF current operator 4 

2 f 

j(r) = v \& r $l(r)<r<*pfo{r) . (A4) 

a, (3=1 ^ 

The perturbing vector potential could in principle induce 
not only a current but also a density modulation. Within 
LRT the induced density can be written in the forrrP^ 

Sn(r,t) = ^2 J dr J d 2 r' XnA r ^ r ' \ T )M^' \t - r) , 

(A5) 



where 



XnAry,t) = -i([n(r,£),j V)])o , (A6) 



with t = {x : y}, is the density-current linear response 
function. For a homogeneous and isotropic system this 
relation takes a much simpler form when written in 
Fourier transform with respect to space and time: 

6n(q, u) = Xnf (<L u)A e (q, u) , (A7) 



with Xn/iQi 00 ) = ((^qi 3- q ))uj- Here we have introduced 
the Kubo productP^ 



((A]B)) u = -i lim 

e— )-0+ 



r+oc 

/ dt 
Jo 



r rf ([i(t),B(o)]) . 

(A8) 

From symmetry arguments Xnj^(Q^) must transform 
as a vector: since q is the only vector available we have 



Xnj*(q,u) = Xnj(q,0j)- 



(A9) 
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where Xnj(q^) dep ends only on the magnitude q of the 
vector q. Thus Eq. (A7) becomes 



8n(q,uj) = Xnj(q^) • (A10) 

This result implies that only longitudinal vector poten- 
tials can produce a density response. 

The evaluation of Xnj(q^) is straightforward. Indeed, 
taking q = qx along the ^-direction, we have 

Xnj(q^) = Xnj*(qX,u) = ((hq',3* q ))u, 

= -{\jq,n-q])o + -Xh(q,u) 



(All) 



where Xnn(q^) and Xl(q,w) are the density- densit}E3 
and longitudinal curr ent- cu r rentP^ response functions, 
respectively. In Eq. (All) we have used the identity 
[A, B] = A 1 "]) 1 ", the following identity valid for Kubo 
products ((i;£))u, = ([A, B]) /oj + i((d t A; and 
Eq. (9) in Ref. [58j Thus, assum ing continuity of the re- 
sponse function, relation ( |A11| ) can be extrapolated to 
the static limit (uo — >• 0) and implies that a static vector 
potential does not give rise to density modulations, since 



lim n Xnj(q^) 

cu— >■() 



lim -Xr, 
uj^o q 



i(q,u) 



(A12) 



For readers who feel uncomfortable with the properties 
of Kubo products we remark that Eq. (All) can also be 



proven explicitly by using the exact eigenstate represen- 
tation for Xnj x (qx : oj) (see Sect. 3.2.3 in Ref. [57]) . 
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